# Replication code for "The Gendered Effects of Violence on Political Engagement"

# Set the appropriate working directory
# setwd("...") 

# Install "foreign" library if you have not already
# install.packages("foreign")

# Load "foreign" library
library(foreign)
# Load "ordinal" library
library(ordinal)

# Import data; ensure that the .dta file is in Stata 12 format (or earlier), otherwise the read.dta() command will not work 
data <- read.dta("Data_Post.dta")

# Create a dummy indicating female respondent status (Woman) 
data$Woman <- ifelse(data$Man == 1, 0, 1)

# Partition the data by experimental condition
data_C_VP <- data[data$Condition != "Identity_Prime",]
data_C_IP <- data[data$Condition != "Violence_Prime",]

# Name experimental conditions
data_C_VP$Treatment <- ifelse(data_C_VP$Condition == "Violence_Prime", 1, 0)
data_C_IP$Treatment <- ifelse(data_C_IP$Condition == "Identity_Prime", 1, 0)

############################################### Main text ###############################################

# Table 1: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_VP)) 
# number of observations
nobs(mod1)
# extract the treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_VP))[[4]][2,] 

# Model 2
# Treatment: Identity Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_IP)) 
# number of observations
nobs(mod2)
# extract the treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_IP))[[4]][2,] 

# Substantive interpretation 
# Constant + Treatment*B1 + Man*B2 + (Treatment x Man)*B3
# Women in control 
-0.095 + 0*-0.183 + 0*0.194 + (0*0)*0.417 # -0.095
# Women in Violence Prime
-0.095 + 1*-0.183 + 0*0.194 + (1*0)*0.417 # -0.278
# Difference
diff(c(-0.095,-0.278)) # -0.183

# Men in control
-0.095 + 0*-0.183 + 1*0.194 + (0*1)*0.417 # 0.099
# Men in Violence Prime
-0.095 + 1*-0.183 + 1*0.194 + (1*1)*0.417 # 0.333
# Difference
diff(c(0.099,0.333)) # 0.234

# Difference between treatment effects
diff(c(-0.183,0.234)) # 0.417

######################################### Online Appendix (OA) #########################################

# OA3: Summary Statistics
# # of respondents
dim(data)[1] # 1125
# # assigned to Control
dim(data[data$Condition=="Control",])[1] # 360
# # assigned to Violence Prime
dim(data[data$Condition=="Violence_Prime",])[1] # 387
# # assigned to Identity Prime
dim(data[data$Condition=="Identity_Prime",])[1] # 378
# % assigned to Control
round(dim(data[data$Condition=="Control",])[1]/1125, 2)*100 # 32
# % assigned to Violence Prime
round(dim(data[data$Condition=="Violence_Prime",])[1]/1125, 2)*100 # 34
# % assigned to Identity Prime
round(dim(data[data$Condition=="Identity_Prime",])[1]/1125, 2)*100 # 34
# # of Bosniaks
dim(data[data$Ethnicity=="Bosniak",])[1] # 608
# # of Croats
dim(data[data$Ethnicity=="Croat",])[1] # 125
# # of Serbs
dim(data[data$Ethnicity=="Serb",])[1] # 392
# % Bosniak
round(dim(data[data$Ethnicity=="Bosniak",])[1]/1125, 2)*100 # 54
# % Croat
round(dim(data[data$Ethnicity=="Croat",])[1]/1125, 2)*100 # 11
# % Serb 
round(dim(data[data$Ethnicity=="Serb",])[1]/1125, 2)*100 # 35
# # of Women
dim(data[data$Man==0,])[1] # 642
# # of Men
dim(data[data$Man==1,])[1] # 483
# % Women
round(dim(data[data$Man==0,])[1]/1125, 2)*100 # 57
# % Men
round(dim(data[data$Man==1,])[1]/1125, 2)*100 # 43
# Mean Age
round(mean(data$Age), 2) # 43.67
# Median Age
median(data$Age) # 43
# Mean Education
round(mean(data$Education), 2) # 4.47
# Median Education
median(data$Education) # 5

# OA5: Balance Tests
# Table OA5.1: Balance Tests, Before Non-Response
# Model 1, Treatment: Violence Prime
summary(mod1 <- glm(Treatment ~ Ethnicity + Man + Age + Education, data=data_C_VP, family=binomial(link="logit")))
nobs(mod1)
# Model 2, Treatment: Identity Prime
summary(mod2 <- glm(Treatment ~ Ethnicity + Man + Age + Education, data=data_C_IP, family=binomial(link="logit")))
nobs(mod2)

# Table OA5.2: Balance Tests, After Non-Response
# Model 1, Treatment: Violence Prime
summary(mod1 <- glm(Treatment ~ Ethnicity + Man + Age + Education, data=data_C_VP[complete.cases(data_C_VP$Engagement_scale),], family=binomial(link="logit")))
nobs(mod1)
# Model 2, Treatment: Identity Prime
summary(mod2 <- glm(Treatment ~ Ethnicity + Man + Age + Education, data=data_C_IP[complete.cases(data_C_IP$Engagement_scale),], family=binomial(link="logit")))
nobs(mod2)

data$Voting_text <- as.character(gsub(1, "Very Unlikely", 
                                      gsub(2, "Fairly Unlikely", 
                                           gsub(3, "Fairly Likely", 
                                                gsub(4, "Very Likely", data$Voting)))))

data$Interest_text <- as.character(gsub(1, "Not At All Interested", 
                                      gsub(2, "A Little Interested", 
                                           gsub(3, "Somewhat Interested", 
                                                gsub(4, "Very Interested", data$Interest)))))

data$Voting_text <- factor(data$Voting_text, levels=c("Very Unlikely","Fairly Unlikely","Fairly Likely", "Very Likely"), ordered=TRUE)
data$Interest_text <- factor(data$Interest_text, levels=c("Not At All Interested","A Little Interested","Somewhat Interested", "Very Interested"), ordered=TRUE)

# OA6: Distribution of Responses Across Experimental Conditions
# Table OA6.1: Distribution of Responses, Voting Question, All Respondents
table(data$Voting_text)
round((table(data$Voting_text)/sum(table(data$Voting_text)))*100)

table(data[data$Condition == "Control", "Voting_text"])
round((table(data[data$Condition == "Control", "Voting_text"])/sum(table(data[data$Condition == "Control", "Voting_text"])))*100)

table(data[data$Condition == "Violence_Prime", "Voting_text"])
round((table(data[data$Condition == "Violence_Prime", "Voting_text"])/sum(table(data[data$Condition == "Violence_Prime", "Voting_text"])))*100)

table(data[data$Condition == "Identity_Prime", "Voting_text"])
round((table(data[data$Condition == "Identity_Prime", "Voting_text"])/sum(table(data[data$Condition == "Identity_Prime", "Voting_text"])))*100)

# Table OA6.2: Distribution of Responses, Voting Question, Male Respondents
dataMen <- data[data$Man == 1,]
table(dataMen$Voting_text)
round((table(dataMen$Voting_text)/sum(table(dataMen$Voting_text)))*100)

table(dataMen[dataMen$Condition == "Control", "Voting_text"])
round((table(dataMen[dataMen$Condition == "Control", "Voting_text"])/sum(table(dataMen[dataMen$Condition == "Control", "Voting_text"])))*100)

table(dataMen[dataMen$Condition == "Violence_Prime", "Voting_text"])
round((table(dataMen[dataMen$Condition == "Violence_Prime", "Voting_text"])/sum(table(dataMen[dataMen$Condition == "Violence_Prime", "Voting_text"])))*100)

table(dataMen[dataMen$Condition == "Identity_Prime", "Voting_text"])
round((table(dataMen[dataMen$Condition == "Identity_Prime", "Voting_text"])/sum(table(dataMen[dataMen$Condition == "Identity_Prime", "Voting_text"])))*100)

# Table OA6.3: Distribution of Responses, Voting Question, Female Respondents
dataWomen <- data[data$Man == 0,]
table(dataWomen$Voting_text)
round((table(dataWomen$Voting_text)/sum(table(dataWomen$Voting_text)))*100)

table(dataWomen[dataWomen$Condition == "Control", "Voting_text"])
round((table(dataWomen[dataWomen$Condition == "Control", "Voting_text"])/sum(table(dataWomen[dataWomen$Condition == "Control", "Voting_text"])))*100)

table(dataWomen[dataWomen$Condition == "Violence_Prime", "Voting_text"])
round((table(dataWomen[dataWomen$Condition == "Violence_Prime", "Voting_text"])/sum(table(dataWomen[dataWomen$Condition == "Violence_Prime", "Voting_text"])))*100)

table(dataWomen[dataWomen$Condition == "Identity_Prime", "Voting_text"])
round((table(dataWomen[dataWomen$Condition == "Identity_Prime", "Voting_text"])/sum(table(dataWomen[dataWomen$Condition == "Identity_Prime", "Voting_text"])))*100)

# Table OA6.4: Distribution of Responses, Interest Question, All Respondents
table(data$Interest_text)
round((table(data$Interest_text)/sum(table(data$Interest_text)))*100)

table(data[data$Condition == "Control", "Interest_text"])
round((table(data[data$Condition == "Control", "Interest_text"])/sum(table(data[data$Condition == "Control", "Interest_text"])))*100)

table(data[data$Condition == "Violence_Prime", "Interest_text"])
round((table(data[data$Condition == "Violence_Prime", "Interest_text"])/sum(table(data[data$Condition == "Violence_Prime", "Interest_text"])))*100)

table(data[data$Condition == "Identity_Prime", "Interest_text"])
round((table(data[data$Condition == "Identity_Prime", "Interest_text"])/sum(table(data[data$Condition == "Identity_Prime", "Interest_text"])))*100)

# Table OA6.5: Distribution of Responses, Interest Question, Male Respondents
dataMen <- data[data$Man == 1,]
table(dataMen$Interest_text)
round((table(dataMen$Interest_text)/sum(table(dataMen$Interest_text)))*100)

table(dataMen[dataMen$Condition == "Control", "Interest_text"])
round((table(dataMen[dataMen$Condition == "Control", "Interest_text"])/sum(table(dataMen[dataMen$Condition == "Control", "Interest_text"])))*100)

table(dataMen[dataMen$Condition == "Violence_Prime", "Interest_text"])
round((table(dataMen[dataMen$Condition == "Violence_Prime", "Interest_text"])/sum(table(dataMen[dataMen$Condition == "Violence_Prime", "Interest_text"])))*100)

table(dataMen[dataMen$Condition == "Identity_Prime", "Interest_text"])
round((table(dataMen[dataMen$Condition == "Identity_Prime", "Interest_text"])/sum(table(dataMen[dataMen$Condition == "Identity_Prime", "Interest_text"])))*100)

# Table OA6.6: Distribution of Responses, Interest Question, Female Respondents
dataWomen <- data[data$Man == 0,]
table(dataWomen$Interest_text)
round((table(dataWomen$Interest_text)/sum(table(dataWomen$Interest_text)))*100)

table(dataWomen[dataWomen$Condition == "Control", "Interest_text"])
round((table(dataWomen[dataWomen$Condition == "Control", "Interest_text"])/sum(table(dataWomen[dataWomen$Condition == "Control", "Interest_text"])))*100)

table(dataWomen[dataWomen$Condition == "Violence_Prime", "Interest_text"])
round((table(dataWomen[dataWomen$Condition == "Violence_Prime", "Interest_text"])/sum(table(dataWomen[dataWomen$Condition == "Violence_Prime", "Interest_text"])))*100)

table(dataWomen[dataWomen$Condition == "Identity_Prime", "Interest_text"])
round((table(dataWomen[dataWomen$Condition == "Identity_Prime", "Interest_text"])/sum(table(dataWomen[dataWomen$Condition == "Identity_Prime", "Interest_text"])))*100)

# OA7: Supplementary Model Results
# Table OA7.1: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, With Pre-Treatment Demographic Controls
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Man + Ethnicity + Age + Education, data=data_C_VP))
nobs(mod1)
# extract the treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman + Ethnicity + Age + Education, data=data_C_VP))[[4]][2,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Man + Ethnicity + Age + Education, data=data_C_IP))
nobs(mod2)
# extract the treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman + Ethnicity + Age + Education, data=data_C_IP))[[4]][2,]

# TABLES OA7.2 AND OA7.3 (MODELS THAT USE SURVEY WEIGHTS) ARE PRESENTED AT THE END OF THIS SCRIPT!

# Table OA7.4: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Voting Question, Without Pre-Treatment Demographic Controls
# Model 1, Treatment: Violence Prime
summary(mod1 <- clm(as.factor(Voting) ~ Treatment*Man, data=data_C_VP, link = "logit"))
# extract the treatment effect for men
summary(clm(as.factor(Voting) ~ Treatment*Woman, data=data_C_VP, link = "logit"))[[5]][4,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- clm(as.factor(Voting) ~ Treatment*Man, data=data_C_IP, link = "logit"))
# extract the treatment effect for men
summary(clm(as.factor(Voting) ~ Treatment*Woman, data=data_C_IP, link = "logit"))[[5]][4,]

# Table OA7.5: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Voting Question, With Pre-Treatment Demographic Controls
# Model 1, Treatment: Violence Prime
summary(mod1 <- clm(as.factor(Voting) ~ Treatment*Man + Ethnicity + Age + Education, data=data_C_VP, link = "logit"))
# extract the treatment effect for men
summary(clm(as.factor(Voting) ~ Treatment*Woman + Ethnicity + Age + Education, data=data_C_VP, link = "logit"))[[5]][4,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- clm(as.factor(Voting) ~ Treatment*Man + Ethnicity + Age + Education, data=data_C_IP, link = "logit"))
# extract the treatment effect for men
summary(clm(as.factor(Voting) ~ Treatment*Woman + Ethnicity + Age + Education, data=data_C_IP, link = "logit"))[[5]][4,]

# Table OA7.6: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Interest Question, Without Pre-Treatment Demographic Controls
# Model 1, Treatment: Violence Prime
summary(mod1 <- clm(as.factor(Interest) ~ Treatment*Man, data=data_C_VP, link = "logit"))
# extract the treatment effect for men
summary(clm(as.factor(Interest) ~ Treatment*Woman, data=data_C_VP, link = "logit"))[[5]][4,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- clm(as.factor(Interest) ~ Treatment*Man, data=data_C_IP, link = "logit"))
# extract the treatment effect for men
summary(clm(as.factor(Interest) ~ Treatment*Woman, data=data_C_IP, link = "logit"))[[5]][4,]

# Table OA7.7: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Interest Question, With Pre-Treatment Demographic Controls
# Model 1, Treatment: Violence Prime
summary(mod1 <- clm(as.factor(Interest) ~ Treatment*Man + Ethnicity + Age + Education, data=data_C_VP, link = "logit"))
# extract the treatment effect for men
summary(clm(as.factor(Interest) ~ Treatment*Woman + Ethnicity + Age + Education, data=data_C_VP, link = "logit"))[[5]][4,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- clm(as.factor(Interest) ~ Treatment*Man + Ethnicity + Age + Education, data=data_C_IP, link = "logit"))
# extract the treatment effect for men
summary(clm(as.factor(Interest) ~ Treatment*Woman + Ethnicity + Age + Education, data=data_C_IP, link = "logit"))[[5]][4,]

# Table OA7.8: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Average Results for All Respondents, Without Pre-Treatment Demographic Controls
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment, data=data_C_VP))
nobs(mod1)
# Model 2, Treatment: Identity Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment, data=data_C_IP))
nobs(mod2)

# Table OA7.9: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Average Results for All Respondents, With Pre-Treatment Demographic Controls
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment + Man + Ethnicity + Age + Education, data=data_C_VP))
nobs(mod1)
# Model 2, Treatment: Identity Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment + Man + Ethnicity + Age + Education, data=data_C_IP))
nobs(mod2)

# OA8: Results Across Respondent Ethnicity
# Table OA8.1: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Bosniak Respondents
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_VP[data_C_VP$Ethnicity == "Bosniak",]))
nobs(mod1)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_VP[data_C_VP$Ethnicity == "Bosniak",]))[[4]][2,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_IP[data_C_IP$Ethnicity == "Bosniak",]))
nobs(mod2)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_IP[data_C_IP$Ethnicity == "Bosniak",]))[[4]][2,]

# Table OA8.2: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Croat Respondents
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_VP[data_C_VP$Ethnicity == "Croat",]))
nobs(mod1)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_VP[data_C_VP$Ethnicity == "Croat",]))[[4]][2,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_IP[data_C_IP$Ethnicity == "Croat",]))
nobs(mod2)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_IP[data_C_IP$Ethnicity == "Croat",]))[[4]][2,]

# Table OA8.3: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Serb Respondents
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_VP[data_C_VP$Ethnicity == "Serb",]))
nobs(mod1)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_VP[data_C_VP$Ethnicity == "Serb",]))[[4]][2,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_IP[data_C_IP$Ethnicity == "Serb",]))
nobs(mod2)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_IP[data_C_IP$Ethnicity == "Serb",]))[[4]][2,]

# Table OA8.4: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Non-Bosniak Respondents
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_VP[data_C_VP$Ethnicity != "Bosniak",]))
nobs(mod1)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_VP[data_C_VP$Ethnicity != "Bosniak",]))[[4]][2,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_IP[data_C_IP$Ethnicity != "Bosniak",]))
nobs(mod2)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_IP[data_C_IP$Ethnicity != "Bosniak",]))[[4]][2,]

# Table OA8.5: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Non-Croat Respondents
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_VP[data_C_VP$Ethnicity != "Croat",]))
nobs(mod1)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_VP[data_C_VP$Ethnicity != "Croat",]))[[4]][2,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_IP[data_C_IP$Ethnicity != "Croat",]))
nobs(mod2)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_IP[data_C_IP$Ethnicity != "Croat",]))[[4]][2,]

# Table OA8.6: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Non-Serb Respondents
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_VP[data_C_VP$Ethnicity != "Serb",]))
nobs(mod1)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_VP[data_C_VP$Ethnicity != "Serb",]))[[4]][2,]
# Model 2, Treatment: Identity Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Man, data=data_C_IP[data_C_IP$Ethnicity != "Serb",]))
nobs(mod2)
# extract treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman, data=data_C_IP[data_C_IP$Ethnicity != "Serb",]))[[4]][2,]

# OA9: Results Across Local Violence Severity
# Table OA9.1: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Controlling for Levels of Violence Severity
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Man + Casualties, data=data_C_VP)) 
nobs(mod1)
# extract the treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman + Casualties, data=data_C_VP))[[4]][2,]
# Model 2, Treatment: Violence Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Man + Log_Casualties, data=data_C_VP)) 
nobs(mod2)
# extract the treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman + Log_Casualties, data=data_C_VP))[[4]][2,] 
# Model 3, Treatment: Violence Prime
summary(mod3 <- lm(Engagement_scale ~ Treatment*Man + Casualty_rate, data=data_C_VP)) 
nobs(mod3)
# extract the treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman + Casualty_rate, data=data_C_VP))[[4]][2,] 
# Model 4, Treatment: Identity Prime
summary(mod4 <- lm(Engagement_scale ~ Treatment*Man + Casualties, data=data_C_IP)) 
nobs(mod4)
# extract the treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman + Casualties, data=data_C_IP))[[4]][2,] 
# Model 5, Treatment: Identity Prime
summary(mod5 <- lm(Engagement_scale ~ Treatment*Man + Log_Casualties, data=data_C_IP)) 
nobs(mod5)
# extract the treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman + Log_Casualties, data=data_C_IP))[[4]][2,] 
# Model 6, Treatment: Identity Prime
summary(mod6 <- lm(Engagement_scale ~ Treatment*Man + Casualty_rate, data=data_C_IP)) 
nobs(mod6)
# extract the treatment effect for men
summary(lm(Engagement_scale ~ Treatment*Woman + Casualty_rate, data=data_C_IP))[[4]][2,] 

# Table OA9.2: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Female Respondents, Across Levels of Violence Severity
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Casualties, data=data_C_VP[data_C_VP$Man == 0,]))
nobs(mod1)
# Model 2, Treatment: Violence Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Log_Casualties, data=data_C_VP[data_C_VP$Man == 0,]))
nobs(mod2)
# Model 3, Treatment: Violence Prime
summary(mod3 <- lm(Engagement_scale ~ Treatment*Casualty_rate, data=data_C_VP[data_C_VP$Man == 0,]))
nobs(mod3)
# Model 4, Treatment: Identity Prime
summary(mod4 <- lm(Engagement_scale ~ Treatment*Casualties, data=data_C_IP[data_C_IP$Man == 0,]))
nobs(mod4)
# Model 5, Treatment: Identity Prime
summary(mod5 <- lm(Engagement_scale ~ Treatment*Log_Casualties, data=data_C_IP[data_C_IP$Man == 0,]))
nobs(mod5)
# Model 6, Treatment: Identity Prime
summary(mod6 <- lm(Engagement_scale ~ Treatment*Casualty_rate, data=data_C_IP[data_C_IP$Man == 0,]))
nobs(mod6)

# Table OA9.3: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Male Respondents, Across Levels of Violence Severity
# Model 1, Treatment: Violence Prime
summary(mod1 <- lm(Engagement_scale ~ Treatment*Casualties, data=data_C_VP[data_C_VP$Man == 1,]))
nobs(mod1)
# Model 2, Treatment: Violence Prime
summary(mod2 <- lm(Engagement_scale ~ Treatment*Log_Casualties, data=data_C_VP[data_C_VP$Man == 1,]))
nobs(mod2)
# Model 3, Treatment: Violence Prime
summary(mod3 <- lm(Engagement_scale ~ Treatment*Casualty_rate, data=data_C_VP[data_C_VP$Man == 1,]))
nobs(mod3)
# Model 4, Treatment: Identity Prime
summary(mod4 <- lm(Engagement_scale ~ Treatment*Casualties, data=data_C_IP[data_C_IP$Man == 1,]))
nobs(mod4)
# Model 5, Treatment: Identity Prime
summary(mod5 <- lm(Engagement_scale ~ Treatment*Log_Casualties, data=data_C_IP[data_C_IP$Man == 1,]))
nobs(mod5)
# Model 6, Treatment: Identity Prime
summary(mod6 <- lm(Engagement_scale ~ Treatment*Casualty_rate, data=data_C_IP[data_C_IP$Man == 1,]))
nobs(mod6)

# CODE FOR TABLES OA7.2 AND OA7.3
rm(list=ls()) 
library(interplot)
library(foreign)
library(ordinal)
library(memisc)
library(nnet)
library(coefplot)
library(effects)
library(pmr)
library(gridExtra)
library(survey)
library(AER)
# Set the appropriate working directory
# setwd("...") 
data <- read.dta("Data_Post.dta")

data$Woman_census <- ifelse(data$Man_census == 1, 0, 1)

data <- data[complete.cases(data$Ethnicity_census),]
data <- data[complete.cases(data$Man_census),]
data <- data[complete.cases(data$Age_census),]
data <- data[complete.cases(data$Education_census),]
data <- data[complete.cases(data$Engagement_scale),]

dim(data) # 1101

data_C_VP <- data[data$Condition != "Identity_Prime",]
data_C_IP <- data[data$Condition != "Violence_Prime",]

data_C_VP$Treatment <- ifelse(data_C_VP$Condition == "Violence_Prime", 1, 0)
data_C_IP$Treatment <- ifelse(data_C_IP$Condition == "Identity_Prime", 1, 0)

# CREATE SURVEY WEIGHTS FOR data_C_VP (CONTROL AND VIOLENCE PRIME)
data_C_VP.svy.unweighted <- svydesign(ids=~1, data=data_C_VP)
# Warning message indicates that the initial survey data are not weighted. We create the survey weights manually using the code below.
gender.dist <- data.frame(Man_census = c(1, 0), Freq = nrow(data_C_VP) * c(0.490566978, 0.509433022))
ethnicity.dist <- data.frame(Ethnicity_census = c("Bosniak", "Croat", "Serb"), Freq = nrow(data_C_VP) * c(0.520299138, 0.160177354, 0.319523508))
age.dist <- data.frame(Age_census = c(1, 2, 3, 4), Freq = nrow(data_C_VP) * c(0.202269683, 0.353700847, 0.26717429, 0.17685518))
education.dist <- data.frame(Education_census = c(1, 2, 3, 4, 5), Freq = nrow(data_C_VP) * c(0.048897384, 0.091729374, 0.214473261, 0.510524395, 0.134375586))
data_C_VP.svy.rake <- rake(design = data_C_VP.svy.unweighted, 
                           sample.margins = list(~Man_census, ~Ethnicity_census, ~Age_census, ~Education_census), 
                           population.margins = list(gender.dist, ethnicity.dist, age.dist, education.dist))
data_C_VP.svy.rake.trim <- trimWeights(data_C_VP.svy.rake, lower=0.3, upper=3, strict=TRUE)
data_C_VP$Weights <- weights(data_C_VP.svy.rake.trim)

# CREATE SURVEY WEIGHTS FOR data_C_IP (CONTROL AND IDENTITY PRIME)
data_C_IP.svy.unweighted <- svydesign(ids=~1, data=data_C_IP)
# Warning message indicates that the initial survey data are not weighted. We create the survey weights manually using the code below.
gender.dist <- data.frame(Man_census = c(1, 0), Freq = nrow(data_C_IP) * c(0.490566978, 0.509433022))
ethnicity.dist <- data.frame(Ethnicity_census = c("Bosniak", "Croat", "Serb"), Freq = nrow(data_C_IP) * c(0.520299138, 0.160177354, 0.319523508))
age.dist <- data.frame(Age_census = c(1, 2, 3, 4), Freq = nrow(data_C_IP) * c(0.202269683, 0.353700847, 0.26717429, 0.17685518))
education.dist <- data.frame(Education_census = c(1, 2, 3, 4, 5), Freq = nrow(data_C_IP) * c(0.048897384, 0.091729374, 0.214473261, 0.510524395, 0.134375586))
data_C_IP.svy.rake <- rake(design = data_C_IP.svy.unweighted, 
                           sample.margins = list(~Man_census, ~Ethnicity_census, ~Age_census, ~Education_census), 
                           population.margins = list(gender.dist, ethnicity.dist, age.dist, education.dist))
data_C_IP.svy.rake.trim <- trimWeights(data_C_IP.svy.rake, lower=0.3, upper=3, strict=TRUE)
data_C_IP$Weights <- weights(data_C_IP.svy.rake.trim)

# CREATE FINAL SURVEY DESIGN OBJECT
design_C_VP <- svydesign(ids=~1, weights=~Weights, data=data_C_VP)
design_C_IP <- svydesign(ids=~1, weights=~Weights, data=data_C_IP)

# Table OA7.2: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Using Survey Weights, Without Pre-Treatment Demographic Controls
# Model 1, Treatment: Violence Prime
summary(OA7.2.1 <- svyglm(Engagement_scale ~ Treatment*Man_census, design_C_VP))
# extract the treatment effect for men
summary(OA7.2.1B <- svyglm(Engagement_scale ~ Treatment*Woman_census, design_C_VP))$coeff[2,]
# extract the number of observations
nobs(OA7.2.1)   
# extract the AIC
OA7.2.1$aic     

# Model 2, Treatment: Identity Prime
summary(OA7.2.2 <- svyglm(Engagement_scale ~ Treatment*Man_census, design_C_IP))
# extract the treatment effect for men
summary(OA7.2.2B <- svyglm(Engagement_scale ~ Treatment*Woman_census, design_C_IP))$coeff[2,]
# extract the number of observations
nobs(OA7.2.2)   
# extract the AIC
OA7.2.2$aic     

# Table OA7.3: The Gendered Effects of Violence and Identity Salience on Desire to Engage in Politics, Using Survey Weights, With Pre-Treatment Demographic Controls
# Model 1, Treatment: Violence Prime
summary(OA7.3.1 <- svyglm(Engagement_scale ~ Treatment*Man_census + Ethnicity_census + Age_census + Education_census, design_C_VP))
# extract the treatment effect for men
summary(OA7.3.1B <- svyglm(Engagement_scale ~ Treatment*Woman_census + Ethnicity_census + Age_census + Education_census, design_C_VP))$coeff[2,]
# extract the number of observations
nobs(OA7.3.1)   
# extract the AIC
OA7.3.1$aic     

# Model 2, Treatment: Identity Prime
summary(OA7.3.2 <- svyglm(Engagement_scale ~ Treatment*Man_census + Ethnicity_census + Age_census + Education_census, design_C_IP))
# extract the treatment effect for men
summary(OA7.3.2B <- svyglm(Engagement_scale ~ Treatment*Woman_census + Ethnicity_census + Age_census + Education_census, design_C_IP))$coeff[2,]
# extract the number of observations
nobs(OA7.3.2)   
# extract the AIC
OA7.3.2$aic 

